ON THE SHOCK WAVE SPECTRUM FOR ISENTROPIC 
GAS DYNAMICS WITH CAPILLARITY 



JEFFREY HUMPHERYS 

Abstract. We consider the stability problem for shock layers in Slem- 
rod's model of an isentropic gas with capillarity. We show that these 
traveling waves are monotone in the weak capillarity case, and become 
highly oscillatory as the capillarity strength increases. Using a spec- 
tral energy estimate we prove that small-amplitude monotone shocks 
are spectrally stable. We also show, through the use of a novel spec- 
tral energy estimate, that monotone shocks have no unstable real spec- 
trum regardless of amplitude; this implies that any instabilities of these 
monotone traveling waves, if they exist, must occur through a Hopf-like 
bifurcation, where one or more conjugate pairs of eigenvalues cross the 
imaginary axis. We then conduct a systematic numerical Evans func- 
tion study, which shows that monotone and mildly oscillatory profiles 
in an adiabatic gas are spectrally stable for moderate values of shock 
and capillarity strengths. In particular, we show that the transition 
from monotone to non-monotone profiles does not appear to trigger any 
instabilities. 



1. Introduction 

We consider Slemrod's model [14^ [36} [37] for a one-dimensional isentropic 
gas with capillarity 



Vt-Uj: = 0, 
(1) , ^/..N _ f!^ 

V / X 



ut+p{v)x = (—] - dv, 

\ V J X 
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where physically, v is the specific volume, u is the velocity in Lagrangian 
coordinates, p{v) is the pressure law for an ideal gas, that is p'{v) < and 
p"{v) > 0, and the coefficient, d > 0, accounting for capillarity strength, is 
constant. This is model is a generalization of the compressible isentropic 
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Navier-Stokes equations, or the j5-system with semi-paraboHc (or real) vis- 
cosity, 



It has recently been shown that viscous shock wave solutions of ^ are spec- 
trally stable for all amplitudes in the case of an adiabatic gas law p{v) = , 
7 G [1,3]; see [H [18]. We remark that this result, together with Mascia & 
Zumbrun's work [311 [30] implies that viscous shocks are asymptotically or- 
bitally stable (hereafter called nonlinearly stable). In this paper, we make 
the first step toward generalizing this work to Slemrod's model by show- 
ing that monotone and mildly oscillatory smooth shock profiles of small to 
moderate amplitude are likewise spectrally stable. 

More generally, we are interested in understanding the degree to which the 
analytic methods used to study shock wave stability in viscous conservation 
laws extend to viscous-dispersive systems. We view Slemrod's model as 
an important test case as it is physically realistic and yet captures some 
of the essential mathematical hurdles found in more extensive models of 
compressible fluid flow. In particular, Slemrod's model is symmetrizable 
and genuinely coupled, having only semi-parabolic diffusion; see [17] for 
details. 

A few notable results in the study of shock wave stability for viscous 
conservation laws include the works of Kawashima |24[ [25l [27] , who proved 
that genuinely coupled symmetrizable systems has a stable essential spec- 
trum, the works of Goodman and others |13[ I32 [ [26 l 120], who proved small- 
amplitude spectral stability for viscous shocks through the use of cleverly 
chosen weighted energy estimates, and the works of Zumbrun and collabora- 
tors j39l [Ml [SH [301 [38] , who proved that spectral stability implies nonlinear 
stability for viscous shocks in conservation laws for both strictly parabolic 
and semi-parabolic viscosities. The missing piece in this overall program is to 
determine whether moderate- and large-amplitude viscous shocks are spec- 
trally stable. Very recently, however, spectral stability for large-amplitude 
shocks for ([2]) was proven in the case of an adiabatic gas [18j, and spectral 
stability was numerically demonstrated for the intermediate range through 
an extensive Evans function study There is some hope that this overall 
strategy will extend to more general systems of viscous conservation laws 
and perhaps even viscous-dispersive models. 

We remark that Kawashima's admissibility results, mentioned above, 
were recently extended to viscous-dispersive (and higher-order) systems [l7j . 
Also, Howard & Zumbrun showed that spectral stability implies nonlinear 
stability for scalar viscous-dispersive conservation laws [161. However, the 
remaining pieces of the general program for viscous-dispersive systems, de- 
scribed above, are still open. 



Vt - Ux 



= 0, 



(2) 




SPECTRUM FOR ISENTROPIC GAS DYNAMICS WITH CAPILLARITY 3 

This paper is organized as follows: In Section O we set the stage by first 
proving the existence of shock profiles for ([T]) through the use of a Lyapunov 
function argument. Then using geometric singular perturbation theory, we 
show that small-amplitude shock profiles converge to the zero-capillarity 
case and are thus monotone. Following that, we use a qualitative ODE 
argument to show that our profiles are monotone for weak capillarity yet 
become highly oscillatory as the capillarity strength d increases. We then 
provide a short estimate on the derivative bounds of the profile, which are 
used later in the stability analysis. Finally, we formulate the integrated 
eigenvalue problem, which makes the stability problem more amenable to 
analysis; see for example 113\ 139] . In Section [3l we generalize the work of 
Matsumura & Nishihara [32j and Barker, Humpherys, Rudd, &: Zumbrun 
[Ij by using a spectral energy estimate to prove that small-amplitude mono- 
tone shocks of ([1]) are spectrally stable. In Section [H we further extend the 
results in ^ and offer a short and novel proof that monotone shocks have no 
unstable real spectrum regardless of amplitude. This restricts the class of 
admissible bifurcations for monotone profiles to those of Hopf-type, where 
one or more conjugate pairs of eigenvalues cross the imaginary axis. The 
approach used here is different than many energy methods in that we use 
a spectral energy estimate that does not appear to have a time-asymptotic 
equivalent, whereas most energy estimates can be performed in either do- 
main. In Section [5l we extend the spectral bounds in [4j to ([1]) by proving 
that high-frequency instabilities cannot occur for adiabatic monotone pro- 
files of any amplitude for d < 1/3. Finally in Section [Gj we carry out a 
systematic numerical Evans function study showing that adiabatic mono- 
tone and mildly oscillatory profiles are spectrally stable for moderate shock 
and capillarity strengths. 

We remark that highly oscillatory profiles in the scalar KDV-Burgers 
model were shown by Pego, Smerka, & Weinstein [3lj to be unstable in cer- 
tain cases. Thus for some, perhaps extreme, parameters, one can reasonably 
expect instabilities to occur in our system as well. It is challenging, how- 
ever, with current technology to explore these extreme cases numerically. 
We plan on exploring this in the future. 

2. Preliminaries 

In this section, we derive the profile ODE and provide a convenient scaling 
for our analysis. We prove the existence of shock profiles for ([T]) through 
the use of a Lyapunov function argument. Then using geometric singular 
perturbation theory, we show that small-amplitude shock profiles converge 
to the zero-capillarity case and are thus monotone. Through a qualitative 
ODE argument, we then show that profiles are monotone for weak capillarity 
yet become oscillatory as the capillarity strength d increases beyond the 
transition point d*. We then provide a short estimate on the derivative 
bounds of the profile, which will be used later in the stability analysis. 
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Finally, we formulate the spectral stability problem and change to integrated 
coordinates making it more amenable to analysis; see for example [131 139j . 

2.1. Shock Profiles. By a shock layer (or shock profile) of ([T]), we mean a 
traveling wave solution 

v{x, t) = v{x — st), 
u{x, t) = u{x — st), 

with asymptotically constant end-states {v,u){±oo) = {v±,u±). Rather by 
translating x ^ x — st, we can instead consider stationary solutions of 

Vt - SVx - Ux = 0, 

Ut - SUx +p{v)rc = ( — ) -dVxxx- 
\ V J X 

Under the rescaling {x,t,u) {—sx,s'^t, —u/s), our system takes the form 

Vt + - = 0, 

^^"^ Ut+Ux + ap{v)a:= (—] - dVxxx-, 

\ V J X 

where a = l/s^. Thus, the shock profiles of ([1]) are solutions of the ordinary 
differential equation 

V — u = U, 

u' + apivy= (^^J -dv'", 

subject to the boundary conditions {v,u){±oo) = {v±,u±). This simplifies 
to 

v' + apiv)' = (^^^ - dv'". 

By integrating from — cxd to x, we get our profile equation, 

v' 

(4) V — + cl{p{v) — p{v-)) = dv" , 

V 

where a is found by setting x = +00, thus yielding the Rankine-Hugoniot 
condition 

(5) a - 



p{v+) -p{v-) 

Without loss of generality, we will assume that < 1;+ < u_ . We remark that 
small-amplitude shocks occur when is close to v- and large-amplitude 
shocks arise to when vj^ nears zero. 

Remark. In the absence of capillarity, that is when d = 0, the profile equa- 
tion (m is of first order, and thus has a monotone solution. As we will 
show, small values of d likewise yield monotone profiles whereas large values 
of d produce oscillatory profiles. We make this precise below. 
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2.2. Adiabatic Gas. Although much of the analysis in this paper holds 
for ideal gases, that is when p'{v) < and p"{v) > 0, our numerical study 
focuses on the special case of an adiabatic gas law, 

(6) p{v) = v~^,j>l, 

together with the rescaling 

(x, t, V, u, a, d) {ex, et,v/e, u/e, ae~'^~^ , e^d), 

where e is chosen so that V- = 1; see [H [18] for more details. This choice 
simplifies our analysis in Section [5] and also gives the Mach number M the 
simplifying form M = Xj ^j^^d. 

2.3. Existence. We prove existence of profiles by the following Lyapunov 
function argument. By writing (jlj) as a first order system, we get 

(7a) v = w, 

(7b) u.' 

where 

(8) (j){v) = v{v — V- + a{p{v) — p{v-)). 

The zero-diffusion case is conservative and has a corresponding Hamiltonian 
that provides us with the needed Lyapunov function. Specifically, let 

(9) E{v,w) = -w^-\ r ^dv. 

2 d V 

Since (j){v) < on then E{v,w) is non-negative for v G [v+,V-]. It 

follows that 

d 

(10) —E{v{x),w{x)) = VE ■ {v', w'f = ^ > 0. 
dx dv 

Hence with diffusion, bounded (homoclinic) orbits at (f+,0) are pulled into 
the minimum (?;_, 0) of E{v, u)) as x — > — oo. Thus there exists a connecting 
orbit from to 



2.4. The Small- Amplitude Limit. We now show that small-amplitude 
shocks of (dJ are monotone and follow the same asymptotic limits as the 
d = case presented in [281 1331 SI HI- We accomplish this by rescaling 
and showing, via geometric singular perturbation theory [IH \T2\ [22] , that 
the profile converges, in the small-amplitude shock limit, to the (monotone) 
non-dispersive case. Thus, monotonicity of small-amplitude shocks of ([1]) is 
implied by the monotonicity of the non-dispersive mentioned above. 

Lemma 2.1. Small- amplitude shocks of ([T]) are monotone for any fixed d. 
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Proof. We scale according to the amplitude e = V- — V-^.. Let v = {v — vo)/e 
and X = ex, where vq = v- — ev-. This frame is chosen so that the end- 
states of the profile are fixed at v± = Tl/2- Additionally, we expand the 
pressure tevm.p{v) and the viscosity term v^^ about V-. Hence (jl]) becomes 



(11) 



e{v - V.) (1 + ap{v.)) + g2 QP>-) (- _ -_)2 ^ ^(gS)^^ _ ^_)3 



By expanding the Rankine-Hugoniot equality, e = a{p{v^) —p{vJ)), about 
we obtain 



(12) 



1 + ap'{v^) 



ap"{v-) 



Substituting (jl2p into (jlip and simplifying gives (recall that v- = 1/2) 



(13) 



ap"{v.\^_^ 



1, 



+ £R{v,v') 



+ e-'dv" 



where R{v,v') = 0(1). Thus, in the e = limit, ([13]) becomes 



(14) 



^, _ ap"{v-)v_ 



which is essentially the same reduction obtained for the viscous Burgers 
equation. Note that the capillarity term vanishes as well, and thus the 
reduction is the same as the zero-capillarity {d = 0) case. 
The slow dynamics of (jlSp take the form 



(15a) 
(15b) 



ew 



w, 

1 

d 



ap"{v- 



\v'-\) + eR{vrv')-^_ 



The fast dynamics, obtained by rescaling x — > x/e, take the form 



(16a) 
(16b) 



w 



ew 

1 

d 



ap"(v_) ,_2 1, w 
I {v^ - ^) + eR{v, v') - — 



We can see from the slow dynamics that solutions will remain on the parabola 
defined by 

_ ap"(y-)v- 2 In 
^ = ^ {v --). 

In addition, we can see from the fast dynamics that any jumps will be 
vertical, that is, v = constant. Since there are no vertical branches, no 
jumps occur and thus it follows that small- amplitude shocks approach the 
solutions for (jl4p . Hence, for sufficiently small amplitudes, the profiles are 
monotone. □ 
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Remark. In the original scale, small- amplitude profiles of ([2]) have the 
asymptotic properties \vx\ = 0{e^) and \vxx\ = \vx\0{e), where e = v- — v+ 
is the amplitude; see \28\ [33] . From the above argument, these asymptotic 
properties hold with our scaling in ([3]) as well. It is also straightforward to 
establish these asymptotic properties directly; see for example Theorem \2.3\ 
below. 



2.5. Classification of Profiles. We show that smooth shock profiles are 
monotone for small values of d and transition to highly oscillatory fronts 
when d gets large; see Figure [J for illustrative examples. The transition 
point between monotone and non-monotone profiles is found to be 

^"^^^ " 4v^_{l+ap'{v_))' 

and in the case of an adiabatic gas with = 1, see Section [221 t^iis becomes 

1 _ M2 
^^^^ * ~ 4(1 - 07) ~ 4(M2-1)' 

where M is the Mach number. In particular as the amplitude approaches 
zero, we have that 1 + ap'(f_) 0, see ()12p . thus making all profiles mono- 
tone regardless of d; this is consistent with the results in Section 12. 4i In the 
large-amplitude limit, we have that a — > and thus — > l/(4u^). Hence, 
for values of d less than l/(4w^), all profiles are monotone, regardless of am- 
plitude, and for d > l/(4f^) a transition from monotone to non-monotone 
occurs for moderate to large amplitude fronts. We have the following: 

Theorem 2.2. Shock profiles of ^ are monotone iffO<d<di,. 

Proof. By a geometric singular perturbation argument very similar to the 
one in Section [2.41 we know that profiles are monotone for sufficiently small 
values of d. When d = d^,, we can show that the local behavior near the 
fixed point (f_, 0) transitions from that of an unstable node to an unstable 
spiral, which is clearly non-monotone. Hence, it suffices to show that the 
profile does not lose monotonicity until d passes through d^: . By linearizing 
((Tj) about we get the system 








(19) y^j = y -j-^+ap'iv-)) 

If monotonicity is lost before d gets to d^ , then for some do < the phase 
curve connects to V- vertically. This would require the vector field near w_ 
to admit a vector in the w direction. However, since 

-w 




dov- 

this cannot happen. Hence the profiles are monotone whenever d < d^. □ 
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0.05 




(c) 

Figure 1. Images of the profiles and their derivatives (left) 
and corresponding phase portraits (right) for an adiabatic 
monatomic gas (7 = 5/3) with = 0.1 and d varying (note 
d:^ « 0.259). We demonstrate (a) monotone profiles with d = 
0.2, (6) non-monotone profiles which are mildly oscillatory 
with d = 2, and (c) non-monotone profiles that are highly 
oscillatory with d = 200. 
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Figure 2. Two phase portraits for an adiabatic monatomic 
gas (7 = 5/3) for = 0.1. The dark hne corresponds to a 
non-monotone profile with d = 5 and the dotted hne is the 
zero-capilarity profile w = (j){v). Note that the dotted curve 
intersects the dark one at its minimum. Hence to bound the 
derivative of the non-monotone profile, we need only bound 
the derivative of the monotone profile. 



2.6. Bounds on \vx\- W now provide bounds on Vx that are used later in 
our analysis. We show that l-O^^I < where e = V- — v+ is the amplitude 
of the profile. This bound holds regardless of capillarity strength, and is 
important for our analysis in Section [H The idea behind the proof follows 
from Figure [21 where we see that the maximum value of \vx\ occurs at the 
point where the profile intersects the zero-capillarity profile. Thus we need 
only find a bound on the zero-capillarity profile. 

Theorem 2.3. Shock profiles of ([T]) satisfy \vx\ < where e = \v^ — v-^-\. 

Proof. Consider the phase portrait of the profile. Let vq denote the point 
that maximizes \vx\- This occurs when w' = in (I7bp . or in other words, 
when w = (j){v), which is the zero-capillarity profile. Hence, the maximum 
point for \vx\ coincides with the zero-capillarity curve, which we can show 
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is bounded above by This follows easily since 

\v— — f+P £^ 

sup|z)^| = sup |</'(f)| < sup \v{v — v^)\ < = — . 

x£S. v€[v+,v-] vG[v+,V-] ^ 4 



□ 



Remark. In the e ^ limit wc ccLTi likewise show that \vxx\ — \'^x 

2.7. Stability problem. By linearizing ([3]) about the profile {v,u), we get 
the eigenvalue problem 

Xv + v' -u' = 0, 



(20) 



Xu + u'- ifiv)vy = ( ^ ) - dv'", 



where f{v) = —ap'{v)—Vx/v'^- We say that a shock profile of ([T]) is spectrally 
stable if the linearized system (j20p has no spectra in the closed deleted 
right half-plane given by P = {3fte(A) > 0} \ {0}, that is, there are no 
growth or oscillatory modes. To show that the essential spectrum is stable, 
we linearize ([3]) about the endstates {v±,u±) and show that the resulting 
constant-coefficient system is stable; see |15] . This was done for general 
viscous-dispersive and higher-order systems in |17j . Thus it suffices to show 
that the point spectrum is also stable. However, since traveling wave profiles 
always have a zero-eigenvalue due to translational invariance, it is often 
difficult to get good uniform bounds in energy estimates. Hence, we use the 
standard technique of transforming into integrated coordinates; see [39t 
S]. This goes as follows: 

Suppose that {v, u) is an eigenfunction of ([20]) with eigenvalue A / 0. 
Then 



u{x) = / u{z)dz, v{x) = / v{z)dz, 

J—OD J—OD 

and their derivatives decay exponentially as x — > oo; see ^9]. Thus, by 
substituting and then integrating, (u, v) satisfies (suppressing the tilde) 

(21a) Xv + v' -u' = 0, 

(21b) Xu + u'- fiv)v' = ^ - dv'" 

V 

This new eigenvalue problem is important because its point spectrum differs 
from that of (j20p only at A = 0. It follows that spectral stability of (j20|) is 
implied by spectral stability of ([2T]) . Hence, we will use (pTj) instead of (pOj) 
in the remainder of our stability analysis. 

3. Small- Amplitude Spectral Stability 

In this section we show that small-amplitude smooth shock profiles are 
spectrally stable. This work generalizes the energy methods in |32l [3] to the 
case of an isentropic gas with capillarity. 
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Theorem 3.1. Small- amplitude shocks of ([T]) are spectrally stable. 

Proof. Suppose that SfteA > 0. Recall that small-amplitude profiles are 
monotone with Vx < and thus also satisfy f{v) > and f'{v) < 0. By 
multiplying ()21bp by the conjugate u/f{v) and integrating in x from — cxd 
to oo, we have 



Xuu 



+ 



u'u 



m 



v'u 



vf{v) 



dv"'u 



Integrating the last three terms by parts and appropriately using (|21ap to 
substitute for u' in the third term gives us 



d 



1 



+ 



f{v) \vf{v) 



1 



u'u + / v{Xv + v') + 



,'|2 



vf{v) 



-v"u + d 



1 



II - 
V u. 



We take the real part and appropriately integrate by parts: 



3?e(A) / 

Jm 



d^e 



+ m 
1 



fiv 



+ 



II -I I 
■V u + 



g{v)\u\'^ + 
1 



u 



l\2 



vf{v) 



where 
(22) 



aiv) 



1 



+ 



1 



II - 

V u 



f{v)J \vf{v) 

Thus, by integrating the last two terms by parts and further simplifying, for 
A > 0, we have 



(23) 



g{v)\u\' + 



< -d^e 



,/|2 



1 



1 



fiv] 



fiv) 



v'u' + 



.'|2 



fiv) 



v'u 



We note that since d > and Vx < 0, then all the terms on the left-hand 
side are non-negative. Moreover, since \vx\ = O(e^) and \vxx\ = |^a;|0(£)> it 
follows that the right-hand side of the above equation is bounded above by 



-2d 



1 



fiv) 



\v'\\u'\ + Cd / e|ua;||f'||n|. 



Thus, by Young's inequality, we have 



9iv)\u\' + 
< -2d 



\u 



l\2 



vf{v 

1 

m 



1 



fii] 



./|2 



4r/i 



+ Vi\V'' 



l\2 



,/|2 



+ c 



£\V. 



,,'|2 



4?72 



+ V2\u\ 
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We can see that for rji > 1 and 7]2, e sufficiently small, the left side dominates 
the right side, which is a contradiction. □ 

4. Monotone large-amplitude shocks 

In this section, we show that monotone profiles have no unstable real 
spectrum. Our proof follows from a novel energy estimate that generalizes 
that of [4J to a general ideal gas law and the addition of a capillarity term. 
This restricts the class of admissible bifurcations for monotone profiles to 
those of Hopf-type, where one or more conjugate pairs of eigenvalues cross 
the imaginary axis. 

Theorem 4.1. Monotone shocks of ([T]) have no unstable real spectrum. 

Proof. Suppose that A € [0,cxd). Since profiles are monotone, we have that 
Vx<^- We multiply ([Mb]) by the conjugate v and integrate in x from — oo 
to oo. This gives 

[ Xuv+ [ u'v - [ f{v)v'v = [ ^-d [ v"'v. 
Jr Jr Jr Jr ^ Jr 

Notice that on the real line, A = A. Thus, we have 

/ Xuv + / u'v — / f{v)v'v = / — — \- d v"v' . 
Jr Jr Jr Jr v 

Using (j21a|) to substitute for Xv in the first term and for u" in the last term, 
we get 

f -,.^[ /■ /• {\v' + v")v , , f 

I u[u —v) + / uv— / j[v)vv= / ha V V. 

Jr Jr Jr Jr ^ J-^ 

Separating terms and simplifying gives 

[ uu' + 2 [ u'v - [ f{v)v'v = X [ ^+ [ ^ + d [ v"v'. 

Jr Jr Jr ^ jr ^ jr 

We further simplify by substituting for u' in the second term and integrating 
the last terms by parts to give, 

[ uu' + 2 [ {Xv + v')v - [ ( f{v) + ^ + ^)v'v+ [ ^-^ = d [ v"v', 
Jr Jr Jr\ v^ V J J^ V J^ 

which yields 

1 uu' + 2X f \v\^ + f (2 + ap'{v) -^v'v+ / ^ = d /" v"v'. 
Jr Jr Jr\ v J v 

By taking the real part (recall that A G [0, cxd)), we arrive at 

This is a contradiction. Thus, there are no positive real eigenvalues for 
monotone shock layers in ([T|) . □ 
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5. High-frequency bounds 



In this section, we prove high-frequency spectral bounds for monotone 
large-amphtude smooth shock profiles. This provides a ceiling as to how 
far along both the imaginary and real axes that one must explore for point 
spectra when doing Evans function computations. Indeed to check for roots 
of the Evans function in the unstable half-plane, say using the argument 
principle, one needs only compute within these bounds. If no roots are 
found therein, then we have a numerical verification of spectral stability. 
We remark that in this section and the next, we depart from the generality 
of an ideal gas, and restrict ourselves to the adiabatic case; see Section \T2[ 
We remark, however, that we could have carried out our analysis for an 
ideal gas as long as V- = 1, which we can achieve by rescaling. We have the 
following lemmata: 



Lemma 5.1. The following identity holds for ei,£2,0 > and ?R.eX > 0; 



(3f?e(A) + |9m(A)|) / + (1 - ei - £3) 



'\2 



U 



(24) 



< 



■uluP + d^ 



1 1 
4^ 2i^ 



,,"|2 



+ I fiv)\vf, 



where C = sxvp\f{v)v\. 



Proof. We multiply (j21b[) by vu and integrate along x from —00 to 00. This 
yields 



A / v\ur + 



vuu + 



.'|2 



f{v)vv'u + d Vxv"u + d vv"u' . 



Taking the real and imaginary parts, adding them together, and noting that 
|JRe(z)| + |3=m(z)| < ^/2\z\, yields 



(5Re(A) + |9m(A)|) / v\u\'^ 

< / v\u\\u'\ + V2 f{ 
Jr. Jr 



2 / '"x\u\ + 



/ki 

Jr 



v]v\v \\u\ 



+ V2d 



\Vx\\v \\u\ + / v\v \\u' 



< El I -uln'P + 



4ei 



v\u\ -|- i 



f{v)\vf + 



29 



+ 



\vx\\u\'^ + 



.."\2 



Vx\\V + £2 V\U\ + 



.'|2 



f{v)v^\u\^ 

262 



V V 



ll\2 



f(v)\v 



l\2 



1 

+ 2 



l^^xl I^P + dF 



1 1 

1^2^2 



II \2 



Rearranging terms yields ([2 



□ 
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Lemma 5.2. The following identity holds for 5ReA > 0; 



(25) 



\u'\^ > 



1 



[fiv)-p'{v)] W 



'|2 



d / w 



Proof. We multiply ()21bp by f' and integrate along x from — oo to oo. This 
yields 

1 



A / uv' 



u'v' 



f{v)\v 



l\2 



-u"v' - d / v"'v' 



R V 



Using (j21a|) on the right-hand side, integrating by parts, and taking the real 
part gives 

../|2 

A / uv + / u'v 



\v'\'^ + ^e{\) [ 
Jr 



+ d \v' 



.11 \2 



In our domain of interest, this yields 



(26) 



A / uv' + I u'v 



1 



/ [f{v)-v'{v)\\v'\'' + d / It;' 



,"|2 



Now we manipulate the left-hand side. Note that 

X [ uv' + [ u'v' = (A + A) [ uv' - [ u{Xv' + v") 
Jr Jr Jr Jr 

= -23f?e(A) [ u'v - [ uu" 
Jr Jr 

= -23f?e(A) [ {Xv + v')v+ [ ju'^. 
Jr Jr 

Hence, by taking the real part we get 



A / uv' + / u'v' 



.'|2 



2Ke(A)^ 



This combines with (j26p to give (j25p . 

Now we prove our high-frequency bounds. 



□ 



Theorem 5.3. For a monotone profile with < d < 1/3, any eigenvalue X 
of (|2ip with nonnegative real part satisfies 

(27) 3f?e(A) + |9m(A)| < 3 + 



5 ' 



where C = sup\f{v)i)\. 

Proof. Combining (j24p and (|25p . we have 

(3f?e(A) + |9m(A)|) / + (1 - ei - £3) 



/(t))|t;r + rf / \v 



."1 2 



< 



v\u\'^ + d^ 



1 1 

4 ^ 2i^ 



/(t))|t; 



/|2 
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Setting 6* = (1 - £1 - £2)72 yields 

(3f?e(A) + |9m(A)|) / ?)|n|2 + (1 - ei - £2)^ 



< 



1 



C 



£1 - £2 

Hence for < d < 1/3, choose £1 = 



{iluP + d^ 



,,"|2 



1 1 

4 ^ 



,,"|2 



1/12 and £2 



1/2 to get ([271). □ 

1, we can show that 



Remark. For an adiabatic gas, p{v) = v~'^ , 7 > 
C < 7; see Thus in the range 7 E [1,3] we can safely bound the un- 
stable spectrum with a half circle of radius 12. This compactifies the region 
of admissible unstable spectrum, thus allowing us to numerically compute 
winding numbers of the Evans function and determine whether shock layers 
are spectrally stable. 

6. Evans function computation 

In this section, we numerically compute the Evans function to determine 
whether any unstable eigenvalues exist in our system. The Evans function 
D(A) is analytic to the right of the essential spectrum and is defined as 
the Wronskian of decaying solutions of (j2ip : see [1]. In a spirit similar to 
the characteristic polynomial, we have that D{X) = if and only if A is 
an eigenvalue of the linearized operator (j2ip . While the Evans function is 
generally too complex to compute analytically, it can readily be computed 
numerically; see [21] and references within. 

Since the Evans function is analytic in the region of interest, we can 
numerically compute its winding number in the right-half plane. This allows 
us to systematically locate roots (and hence unstable eigenvalues) within. As 
a result, spectral stability can be determined, and in the case of instability, 
one can produce bifurcation diagrams to illustrate and observe its onset. 
This approach was first used by Evans and Feroe [lO] and has been applied 
to various systems since; see for example [341 [2| ISllG]. 

6.1. Numerical Setup. We begin by writing 
W = A{x, X)W, where 



211) as a first-order system 



(28) A{x,X) 



1 


A 


1 


^ 








1 














1 



w 



\X/d X/d h/d -{dvy^J 



V'J 



X/v. Note that eigenvalues of ^1 



and h = h{v, A) := 1 -|- ap'{v) + Vx/v^ 
correspond to nontrivial solutions of PF(x) for which the boundary condi- 
tions VF(iboo) = are satisfied. We remark that since v is asymptotically 
constant in then so is A(x,A). Thus at each end-state, we have the 
constant-coefficient system 

(29) W' = A^{X)W, A±(A) := lim A[x,X). 

x— >±oo 
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Hence solutions that satisfy the needed boundary condition must emerge 
from the 2-dimensional unstable manifold Wi{x) A VFg"!^) at x = — oo and 
also the 2-dimensional stable manifold W^{x) A W^{x) at = cxd. In other 
words, eigenvalues of (j2ip correspond to the values of A for which these two 
manifolds intersect, or more precisely, when -D(A) = 0, where 

D{X) := {W,~ A A A W+)\,^o = detiW^W^W+Wt)i,^o- 

We cannot naively produce the stable and unstable manifolds numeri- 
cally. Indeed with two exponential growth and decay modes, problems with 
stiffness arise. Hence, we use the compound-matrix method to analytically 
track the stable and unstable manifolds; see [31 IH El El [21]. Specifically we 
lift A{x,X) into the exterior-product space A^(C^) « to get 








1 





-1 


















1 


A 












X/d 


h/d 


-idv)-^ 





A 


1 


















1 









-x/d 








h/d 


-{dvy^ 


1 




v 





-X/d 





-X/d 





-{dv) 





We then consider single trajectories W±{x) of the "lifted" problem 

W' = A^^\x,X)W 

on each side corresponding to the simple dominant growth and decay modes 
at the left and right end states, respectively. These trajectories correspond 
to the 2-forms W^{x) A W2{x) and Wj^{x) A W^{x), and can be effectively 
wedged together when they meet at zero; see [^ for an excellent overview 
of this method. 

As an alternative, we consider the adjoint formulation of the Evans func- 
tion [35l[5]. Specifically, we integrate the trajectory W+ along the largest 
growth mode of the adjoint ODE 

(30) w' = -A^^\x,xyw. 

starting at x = oo. We then define the (adjoint) Evans function to be 
D+{X) := {W+ ■ W-)\x=o- Note that W-^- corresponds to the orthogonal 
complement of the 2-form W^{x) A W^{x) and so orthogonality of W+ and 
W- corresponds to intersection of the stable and unstable manifolds. 

To further improve the numerical efficiency and accuracy of the shooting 
scheme, we rescale W and W to remove exponential growth/decay at infinity, 
and thus eliminate potential problems with stiffness. Specifically, we let 
W{x) = e'^ ^V{x), where fi~ is the largest growth rate of the unstable 
manifold at x = — oo, and we solve instead V'{x) = {A^'^^x, A) — fi~I)V{x). 

(2) 

We initialize V{x) at x = — oo as eigenvector r~ of A_ [X) corresponding 

to fi~ . Similarly, it is straightforward to rescale and initialize W{x) at 
X = oo. This method is known to have excellent accuracy [71 [El [H [6l [2T1 14]: 
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Figure 3. Evans function output for semi-circular contour 
of radius 12 with d = 0.45 and (a) v+ = 0.65, (b) v+ = 0.45, 
(c) v+ = 0.35, (d) v+ = 0.25, (e) v+ = 0.20, and (f) v+ = 
0.15. Although the contours wrap around the origin as the 
shock strength increases, they clearly have winding number 
zero, thus demonstrating spectral stability. 



in addition, the adaptive refinement gives automatic error control. Finally, 
in order to maintain analyticity, the initial eigenvectors r~{\) are chosen 
analytically using Kato's method; see [23l pg. 99] and also [9l [6t [T9]. 
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0.08 r 




Figure 4. Evans function output of a semi-circular contour 
with d = 0.75 and v+ G [0.20,0.80]. As the shock strength 
increases, the contours get closer to the origin and begin to 
wrap around it. In the small shock limit, the contour drifts 
away from the origin and gets smaller. 

6.2. Numerical Experiments. We truncate the domain to a sufficiently 

large interval [L_ , Lj^] in order to do numerical computation. Some care 
needs to be taken, however, to make sure that we go out far enough to 
produce good results. Our experiments, described below, were primarily 
conducted using L± = ±25, but for weaker shocks we had to go out as far 
as L± = ±50. For highly oscillatory profiles, very large values of L± are 
needed because the (under-damped) decay rate can be small. To compute 
the profile, we used Matlab's bvp4c routine, which is an adaptive Lobatto 
quadrature scheme. 

Our experiments were carried out uniformly on the range 

{v+,d) £ [0.10,0.80] X [0.05,0.80], 

with 7 = 1.4 (diatomic gas). In terms of Mach number, this corresponds 
roughly to 1.15 < M < 5, which covers the supersonic range and goes into 
the hypersonic regime. Indeed M ~ 5 may even go beyond the physical range 
of the model. For each {v-^,d) on our grid, we computed the Evans function 
along a semi-circular contour in the right-half plane of radius 12 centered 
at the origin. Recall that for d < 1/3, this contains the admissible region of 
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0.2r 




Figure 5. Evans function output of a semi-circular contour 
with = 0.25 and d G [0.15,0.80]. As the d decreases, the 
contours get larger and more spread out. 

unstable spectrum from our high-frequency bounds. The ODE calculations 
for individual values of A were carried out using Matlab's ode45 routine, 
which is the adaptive 4th-order Runge-Kutta-Fehlberg method (RKF45). 
Typical runs involved between 100 and 700 mesh points, with error tolerance 
set to AbsTol = le-6 and RelTol = le-8. Values of A were varied on the 
semi-circular contour with 70 points in the first quadrant, 40 on the arc and 
30 along the imaginary axis, and then reflected along the real axis due to 
the conjugate symmetry of the Evans function, that is, -D(A) = -D(A). 

In Figure [3l we see a typical run for increasing Notice that the contour 
wraps around the origin as the shock strength increases. Thus it is difficult 
to conclude stability in the strong shock limit; this is a topic for future 
consideration. Note that the graph gets farther away from the origin as 
the shock strength decreases, thus strongly suggesting stability in the small- 
amplitude limit. In Figure [H we see this effect more clearly. In Figure [5l 
we hold the shock strength fixed and vary d. As d approaches zero, we see 
the contour getting larger and more spread out. Otherwise output does not 
seem to vary much in d, at least in our region of interest. 

The actual parameter values computed were 

iv+,d) G {0.10, 0.15, . . . , 0.80} x {0.05, 0.10, . . . , 0.80}; 
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0.1 - 



gl , 1 1 1 1 , , , 1 1 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

V 

+ 

Figure 6. Dots correspond to runs with parameters {v+,d). 
The upward increasing curve corresponds to the critical value 
d^: between monotone and oscillatory shock profiles. 



see Figure [6l In total 240 runs were conducted, all of which had winding 
number zero. This effectively demonstrates spectral stability for monotone 
and nearly monotone profiles with d < 1/3 and strongly suggests spectral 
stability elsewhere in our region of study. Indeed the output is strikingly 
similar throughout. Nonetheless, for d >> 1 our profile becomes highly 
oscillatory and so it is not unreasonable to expect an instability to occur in 
the extreme. This is a good direction for future work. 

7. Discussion and Open Problems 

We note that (1280 blows up as f + — > 0, and moreover the eigenvalues get 
far apart, thus causing extreme stiffness. Hence we have numerical difficul- 
ties for strong shocks, e.g, M >> 5. Difficulties also arise for both large and 
small values of d. In particular the profile becomes highly oscillatory and 
numerically intractable for very large values of d, and as d — > we likewise 
have that ()28p blows up. Nonetheless, we may be able to demonstrate sta- 
bility as d ^ analytically as a singular limit of the d = case, which is 
stable; see [18] . 

Slemrod's model is an ideal system for further investigation. Not only 
is it physically relevant, and in some sense a canonical viscous-dispersive 
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system, but it also pushes the boundaries of current numerical methods. 
While this model has nice features such as monotone profiles, it also has 
highly complex and numerically taxing obstacles such as highly oscillatory 
profiles and large spectral separation between modes of (|28p in the extreme 
parameter regime. We intend to study this model further. 
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